function ay = heightBasedControl(hf, h, tgo, D, speed, fpa, ay0)
    global kh
    dhc = - (D/m) * 2 * 7110 / V;
    dh = speed * sin(fpa);
	aQEGC = kh * (dhc - dh);
    hf_est = h + dh * tgo;
    if tgo > 800
		ah = 0.2 / tgo / tgo * (hf - hf_est);
	elseif tgo > 600
		ah = (0.2 + 3.8 * (tgo - 600) / 200) / tgo / tgo * (hf - hf_est);
	elseif tgo > 10
		ah = 4.0 / tgo / tgo * (hf - hf_est);
	elseif tgo > 1
		ah = (0.05 + 3.95 * (tgo) / 9) / tgo / tgo * (hf - hf_est);
	else
		ah = 0.05 * tgo * (hf - hf_est);
    end
    ay = ay0 + aQEGC + ah;
end